home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.4 / xfm / fm.c next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  6.9 KB  |  327 lines
  1. /* 
  2.  * fm.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * F(ast)M(oments).C
  11.  *
  12.  * routines to extract, from the input image, values of the coefficients 
  13.  * for the Hatamian filter array and to evaluate various moments and
  14.  * related quantities
  15.  */
  16.  
  17. #include "xfm.h"
  18.  
  19. #define    ZERO        0
  20. #define    NOLOOP        ZERO
  21.  
  22. #undef    DUMP_IMG
  23. #undef    F_DEBUG
  24. #undef    AP_DEBUG
  25. #define    DEBUG
  26.  
  27. /*
  28.  * filter_recursion()
  29.  *   DESCRIPTION:
  30.  *     extract filter array from image
  31.  *   ARGUMENTS:
  32.  *     lx: leftmost x (int)
  33.  *     rx: rightmost x (int)
  34.  *     nr: number of rows (int)
  35.  *     nc: number of columns (int)
  36.  *     y: output array for filter (float **)
  37.  *     xdim: x size of moment matrix (int)
  38.  *     ydim: y size of moment matrix (int)
  39.  *     imgIn: image (Image *)
  40.  *   RETURN VALUE:
  41.  *     number of mode parameters
  42.  */
  43.  
  44. void
  45. filter_recursion (lx, rx, nr, nc, y, xdim, ydim, imgIn)
  46.      int lx, rx;
  47.      int nr, nc;
  48.      float *y[];
  49.      int xdim, ydim;
  50.      Image *imgIn;
  51. {
  52.   int p, q;
  53.   int i, j;
  54.   int imij;
  55.   double hz1, hz2, hz3, hz4;
  56.  
  57.  
  58.   printf ("\n...executing recursion to obtain Y - matrix...\n");
  59.  
  60.   printf ("\n...row:");
  61.   for (i = 1; i < nr; i++) {
  62.     if (i % 50 == 0)
  63.       printf ("...%3d", i);
  64.  
  65.     hz1 = hz2 = hz3 = hz4 = 0.0;
  66.     for (j = 0; j < nc; j++) {
  67.       if ((imij = (int) getpixel (j, i, imgIn)) != 0)
  68.         imij /= imij;
  69.       //if( (imij = inq_index(i, j)) != 0 )   imij /= imij;
  70.  
  71.       hz4 += hz3;
  72.       hz3 += hz2;
  73.       hz2 += hz1;
  74.       hz1 += (double) imij;
  75.     }
  76.  
  77.     for (p = ydim - 1; p >= 1; p--) {
  78.       for (q = 0; q < xdim; q++)
  79.         *(y[q] + p) += *(y[q] + (p - 1));
  80.     }
  81.     *(y[0] + 0) += (float) hz1;
  82.     *(y[1] + 0) += (float) hz2;
  83.     *(y[2] + 0) += (float) hz3;
  84.     *(y[3] + 0) += (float) hz4;
  85.   }
  86.  
  87. #ifdef F_DEBUG
  88.   printf ("\n...filter output (Y-matrix):\n");
  89.   for (p = 0; p < ydim; p++) {
  90.     for (q = 0; q < xdim; q++) {
  91.       printf ("  y[%d][%d] = %f ", p, q, *(y[p] + q));
  92.       printf ("\n");
  93.     }
  94.   }
  95. #endif
  96. }
  97.  
  98.  
  99. #if 0
  100. /*
  101.  * AP code to evaluate moments, given filter output in 2d array Y
  102.  *
  103.  * implementation of M. Hatamian's code involving matrix 
  104.  * element operations in form of matrix multiplication:
  105.  *
  106.  *   m = cm*y*cm';   CM': transposed CM, *: matrix mult.
  107.  */
  108. void
  109. ap_eval_mom (y, m, mdim)
  110.      float *y[], *m[];
  111.      long mdim;
  112. {
  113.   int i, j;
  114.   float *y_out;
  115.  
  116.   long mmdd;
  117.   long ap_dim;
  118.   long ap_y, ap_cm, ap_m;
  119.   long handle_1;
  120.  
  121.   mmdd = mdim * mdim;
  122.  
  123. #ifdef AP_DEBUG
  124.   printf ("\n...input data:\n");
  125.   for (i = 0; i < (int) mdim; i++) {
  126.     for (j = 0; j < mdim; j++) {
  127.       printf (" %07.4f ", y[i][j]);
  128.       if ((j + 1) % 4 == 0)
  129.         printf ("\n");
  130.     }
  131.     printf ("\n");
  132.   }
  133. #endif
  134.  
  135. /*
  136.  *  allocate memory
  137.  */
  138.   if ((y_out = (float *) calloc (mmdd, sizeof (float))) == NULL)
  139.       fail_alloc ("y", 1);
  140. /*
  141.  *  allocate AP memory
  142.  */
  143.   ap_dim = (long) ipaloc (1L);
  144.   ap_y = (long) ipaloc (mmdd);
  145.   ap_m = (long) ipaloc (mmdd);
  146.   ap_cm = (long) ipaloc (mmdd);
  147.  
  148. #ifdef AP_DEBUG
  149.   apmdia ();
  150. #endif
  151.  
  152.   printf ("\n...transfer data to array processor\n");
  153. /*
  154.  * load data onto AP
  155.  */
  156.   apputa ((long) y, ap_y, mmdd);
  157.   apputa ((long) cm, ap_cm, mmdd);
  158.   apputv ((float) mdim, ap_dim);
  159.  
  160.  
  161. #ifdef AP_DEBUG
  162. /*
  163.  * check whether appropriate values have been transferred
  164.  */
  165.   printf ("\n...input matrix, as read back:\n");
  166.   apgeta (ap_y, (long) y_out, mmdd);
  167.   for (i = 0; i < (int) mmdd; i++) {
  168.     printf ("%07.4f  ", *(y_out + i));
  169.     if ((i + 1) % 4 == 0)
  170.       printf ("\n");
  171.   }
  172.   apputa ((long) y_out, ap_y, mmdd);
  173.  
  174. #endif
  175.  
  176. /*
  177.  * instructions to execute the operations of interest
  178.  */
  179.   apmac (&handle_1);
  180.   rmmul (ap_cm, ap_dim, ap_y, ap_dim, ap_y, ap_dim,
  181.          ap_dim, ap_dim, ap_dim);
  182.   rtrn2 (ap_cm, ap_dim, ap_dim);
  183.   rmmul (ap_y, ap_dim, ap_cm, ap_dim, ap_m, ap_dim,
  184.          ap_dim, ap_dim, ap_dim);
  185.   apendm (handle_1, NOLOOP);
  186.  
  187. /* pass handles to routines  */
  188.   apgo (handle_1);
  189.   printf ("\n\n...ckeck DONE flag (do not wait if not set): ");
  190.   printf (" %d\n", ipdtst ());
  191.  
  192. /* output is now available and may be fetched */
  193.   apgeta (ap_m, (long) m, mmdd);
  194.   printf ("\n...moment output:\n");
  195.   for (j = 0; j < (int) mdim; j++) {
  196.     for (i = 0; i < (int) mdim; i++) {
  197.       printf ("  %07.4f", m[i][j]);
  198.       if ((i + 1) % (int) mdim == 0)
  199.         printf ("\n");
  200.     }
  201.   }
  202.   free (y_out);
  203. }
  204.  
  205. #endif
  206.  
  207. /*
  208.  * eval_mom()
  209.  *   DESCRIPTION:
  210.  *     evaluate moments, given filter output in 2d array Y
  211.  *     implementation of M. Hatamian's code involving matrix
  212.  *     element operations in form of matrix multiplication:
  213.  *     m = cm*y*cm';    CM': transposed CM, *: matrix mult.
  214.  *   ARGUMENTS:
  215.  *     y: input matrix (float **)
  216.  *     m: output matrix (float **)
  217.  *     mdim: number of rows, cols in matrices (int)
  218.  *   RETURN VALUE:
  219.  *     none
  220.  */
  221.  
  222. void
  223. eval_mom (y, m, mdim)
  224.      float *y[], *m[];
  225.      int mdim;
  226. {
  227.   int i, j, k;
  228.  
  229. #ifdef DEBUG
  230.   printf ("\n...coefficient arrays\n");
  231.   printf ("...cm:\n");
  232.   for (j = 0; j < mdim; j++) {
  233.     for (i = 0; i < mdim; i++) {
  234.       printf (" %4.2f ", cm[j][i]);
  235.     }
  236.     printf ("\n");
  237.   }
  238.   printf ("...cmt:\n");
  239.   for (j = 0; j < mdim; j++) {
  240.     for (i = 0; i < mdim; i++) {
  241.       printf (" %4.2f ", cmt[j][i]);
  242.     }
  243.     printf ("\n");
  244.   }
  245. #endif
  246.  
  247.  
  248. /* matrix multiplications */
  249.   for (j = 0; j < mdim; j++) {
  250.     for (i = 0; i < mdim; i++) {
  251.       *(m[j] + i) = (float) 0.0;
  252.       for (k = 0; k < mdim; k++) {
  253.         *(m[j] + i) += cm[j][k] * (*(y[k] + i));
  254.       }
  255.     }
  256.   }
  257.  
  258.   for (j = 0; j < mdim; j++) {
  259.     for (i = 0; i < mdim; i++) {
  260.       *(y[j] + i) = (float) 0.0;
  261.       for (k = 0; k < mdim; k++) {
  262.         *(y[j] + i) += *(m[j] + k) * cmt[k][i];
  263.       }
  264.     }
  265.   }
  266.  
  267.  
  268.   for (j = 0; j < mdim; j++) {
  269.     for (i = 0; i < mdim; i++) {
  270.       *(m[j] + i) = *(y[j] + i);
  271.     }
  272.   }
  273. }
  274.  
  275.  
  276.  
  277.  
  278. /*
  279.  * evaluate central moments
  280.  */
  281. /*
  282.  * central_moments()
  283.  *   DESCRIPTION:
  284.  *     evaluate central moments
  285.  *   ARGUMENTS:
  286.  *     m: input matrix (float **)
  287.  *     mu00: central moment (float *)
  288.  *     mu11: central moment (float *)
  289.  *     mu02: central moment (float *)
  290.  *     mu20: central moment (float *)
  291.  *     xc: centroid (float *)
  292.  *     yc: centroid (float *)
  293.  *     r_gyr: radius of gyration (float *)
  294.  *     dent: eccentricity (float *)
  295.  *   RETURN VALUE:
  296.  *     none
  297.  */
  298.  
  299. void
  300. central_moments (m, mu00, mu11, mu02, mu20, xc, yc, r_gyr, dent)
  301.      float **m;
  302.      float *mu00, *mu11, *mu02, *mu20;
  303.      float *xc, *yc, *r_gyr, *dent;
  304. {
  305.   float vmu00;
  306.  
  307. /*
  308.  * centroid 
  309.  */
  310.   *xc = (*(m[1] + 0) / *(m[0] + 0));
  311.   *yc = (*(m[0] + 1) / *(m[0] + 0));
  312.   printf ("\ncentroid:\n");
  313.   printf ("...x_bar = %f  y_bar = %f\n", *xc, *yc);
  314.  
  315. /*
  316.  * central moments and radius of gyration
  317.  */
  318.   printf ("...m00 = %f\n", *(m[0] + 0));
  319.   vmu00 = *mu00 = *(m[0] + 0);
  320.   *mu11 = (*(m[1] + 1) - vmu00 * (*(m[1] + 0) / vmu00) * (*(m[0] + 1) / vmu00));
  321.   *mu20 = (*(m[2] + 0) - vmu00 * (*(m[1] + 0) / vmu00) * (*(m[1] + 0) / vmu00));
  322.   *mu02 = (*(m[0] + 2) - vmu00 * (*(m[0] + 1) / vmu00) * (*(m[0] + 1) / vmu00));
  323.  
  324.   *r_gyr = (float) sqrt ((*mu20 + *mu02) / (vmu00));
  325.   *dent = ((*mu02 - *mu20) / (*mu20 + *mu02));
  326. }
  327.